clc;
clear;

%% main主脚本
tic;
startTime = cputime;
pde = process_pde_NS_s(1);
vector_h = [1/4 1/8 1/16 1/32];
basis_type1 = "P2";
basis_type2 = "P1";
iter_L = 8;

u_Linf = zeros(size(vector_h,2),1);
u_L2 = zeros(size(vector_h,2),1);
u_H1 = zeros(size(vector_h,2),1);
p_Linf = zeros(size(vector_h,2),1);
p_L2 = zeros(size(vector_h,2),1);
p_H1 = zeros(size(vector_h,2),1);

for i = 1:size(vector_h,2)
    h = vector_h(i);
    mesh = process_mesh(pde, h);
    fs1 = process_fespace(mesh, basis_type1);
    fs2 = process_fespace(mesh, basis_type2);

    % process_showmesh(mesh, fs1);
    % process_showmesh(mesh, fs2);

    A12 = pde.nu*term_grad_U_grad_V(mesh, fs1);
    A13 = term_P_div_V(mesh, fs1, fs2);
    O1 = sparse(fs2.Nb, fs2.Nb);
    O2 = sparse(fs2.Nb, 1);
    A = [A12, -A13;
        -A13', O1];
    b1 = term_F_V(mesh, fs1, pde.f1, pde.f2);
    Un = sparse(2*fs1.Nb, 1);

    for l = 1:iter_L
        A111 = term_Uold_grad_U_V(mesh, fs1, Un);
        A112 = term_U_grad_Uold_V(mesh, fs1, Un);
        AN = blkdiag(A111+A112, O1);
        Al = A + AN;
        bN = term_Uold_grad_Uold_V(mesh, fs1, Un);
        bl = [b1 + bN; O2];
        [Al, bl] = process_boundary_NS_s(Al, bl, mesh, fs1, fs2, pde);
        SOL = Al\bl;
        Un = SOL(1:2*fs1.Nb);
    end

    % process_field(fs1, fs2, SOL);

    [u_Linf(i,1), p_Linf(i,1)] = process_error_NS_s(mesh, fs1, fs2, pde, SOL, "Linf");
    [u_L2(i,1), p_L2(i,1)] = process_error_NS_s(mesh, fs1, fs2, pde, SOL, "L2");
    [u_H1(i,1), p_H1(i,1)] = process_error_NS_s(mesh, fs1, fs2, pde, SOL, "H1");

    fprintf("已完成\t%d/%d\n", i, size(vector_h,2));
end

process_rate(u_Linf, vector_h);
process_rate(u_L2, vector_h);
process_rate(u_H1, vector_h);
process_rate(p_Linf, vector_h);
process_rate(p_L2, vector_h);
process_rate(p_H1, vector_h);

TOCtime = toc;
CPUtime = cputime - startTime;
disp(["TOC time(s)", num2str(TOCtime)]);
disp(["CPU time(s)", num2str(CPUtime)]);

